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A local dynamic model for large eddy simulation 

By S. Ghosal, T. S. Lund AND P. Moin 
1. Motivation and objectives 

The dynamic model (Germano et al 1991) is a method for computing the co- 
efficient C in Smagorinsky’s (1963) model for the subgrid-scale stress tensor as a 
function of position from the information already contained in the resolved velocity 
field rather than treating it as an adjustable parameter. There are two advantages 
to this. Firstly, it gives a systematic procedure for computing a flow about which 
we have no prior experience and, therefore, cannot properly adjust the parameter 
C. Secondly, in an inhomogeneous flow, the optimum choice of C may be different 
at different points in the flow, and one does not expect the entire flow to be rep- 
resented by a single constant. The same consideration applies to flows undergoing 
a transition to turbulence or, more generally, to flows whose statistical properties 
are changing with time. In the traditional approach, one needs to introduce further 
ad hoc assumptions, such as wall damping functions or a prescription to reset the 
value of C from zero to a finite number as the flow undergoes a transition to turbu- 
lence. In contrast, inhomogeneous and statistically unsteady flows can be handled 
very naturally in the context of the dynamic model since C is now a function of 
position and time. Though the dynamic model lacked the full generality necessary 
to handle general turbulent flows with no homogeneous directions, the method had 
some important successes. 

The basic formalism behind the method is summarized below following Germano 
et al, (1991). The equations of LES can be thought of as a filtered form of the 
Navier-Stokes equations where the filtering serves to remove fluctuations on length- 
scales smaller than the computational grid. The effect of the unresolved eddies 
on the large eddies is then manifested through the Reynolds stress term r tJ = 

UiUj — UiUj where the bar denotes some grid-level filtering operation on a given 
field vK x ); 


t/>(x 


>=/ 


G 0 (if;y)4’(y)dy. 


The filtering kernel Go(x,y) has a ‘filter width’ equal to the grid spacing A of the 
LES. To compute C one first introduces a ‘test’ filtering operation on the large-eddy 
field that is denoted by the symbol ‘ 


0(x) = J G(x,y)V>(y)dy, 


where G(x,y) is any kernel that serves to damp all spatial fluctuations shorter than 
some characteristic length A > A and x, y are position vectors. The equations for 
the test-filtered field contain the Reynolds stress term Tij = uju] — u,Uj. Both T,j 
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and Tij are unknown in LES; however, the two are related by the identity (Germano 
1992) 

Lij = Tij - t^. ( 1 ) 

Here the Leonard term L i} = u t u 3 - %Hj is computable from the large-eddy field. 
Finally, it is assumed that a seeding law is operative and, therefore, the Reynolds 
stress at the grid and test levels may be written as 

Tij - ±T kk 6ij = — 2CA 2 |3|3 ii . (2) 

a nd ^ ^ 

Tij - \T kk Sij = -2CA 2 |f|f,„ (3) 

respectively. The model coefficient ‘C’ in (2) and (3) need not be the same. The 
prescription for determining C described below can be generalized to obtain both 
coefficients (Moin 1991). In what follows, ‘C” is taken to be the same in (2) and (3) 
for simplicity. On substituting (2) and (3) in (1), an equation for determining C is 
obtained: 

Lij - ^Lkkhj = oiijC - PijC (4) 

where ^ ^ ^ 

oiij = — 2A 2 |S|Sij, 

ftij = -2A 2 |S|^ t> . 

Since C appears inside the filtering operation, equation (4) is a system of five 
independent integral equations involving only one function C. In previous formu- 
lations (Germano ei al 1991, Moin et ai 1991, Lilly 1992), one simply ignored the 
fact that C is a function of position and took C out of the filtering operation as if it 
were a constant. This ad hoc procedure cannot even be justified a posteriori because 
the C field computed using this procedure is found to be a rapidly varying function 
of position (Moin 1991). One of the objectives of this research is to eliminate this 
mathematical inconsistency. 

The C obtained from equation (4) can be either positive or negative. A negative 
value of C implies a locally negative eddy- viscosity, which in turn implies a flow of 
energy from the small scales to the resolved scales or back-scatter. It is known from 
direct numerical simulation (DNS) data (Piomelli et al. 1991) that the forward 
and reverse cascade of energy in a turbulent flow are typically of the same order 
of magnitude with a slight excess of the former accounting for the overall trans- 
fer of energy from large to small scales. The presence of back-scatter, therefore, 
is a desirable feature of a subgrid-scale model. However, when the C computed 
from (4) is used in a large-eddy simulation, the computation is found to become 
unstable. The instability can be traced to the fact that C has a large correlation 
time. Therefore, once it becomes negative in some region, it remains negative for 
excessively long periods of time during which the exponential growth of the local 
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velocity fields, associated with negative eddy- viscosity, causes a divergence of the 
total energy. Though this issue of stability remained unresolved, a way around the 
problem was found if the flow possessed at least one homogeneous direction. Pre- 
vious authors (Germano et ai 1991, Moin et al. 1991, Cabot and Moin 1991) have 
used an ad hoc averaging prescription to stabilize the model. The disadvantages of 
this are: (a) It is based on an ad hoc procedure, (b) The prescription can only be 
applied to flows that have at least one homogeneous direction, thus excluding the 
more challenging flows of engineering interest, (c) The prescription for stabilizing 
the model makes it unable to represent back-scatter. The present research attempts 
to eliminate these deficiencies. 

2. Accomplishments 

In the next section, a variational formulation of the dynamic model is described 
that removes the inconsistency associated with taking C out of the filtering opera- 
tion. This model, however, is still unstable due to the negative eddy- viscosity. Next, 
three models are presented that are mathematically consistent as well as numeri- 
cally stable. The first two are applicable to homogeneous flows and flows with at 
least one homogeneous direction, respectively, and are, in fact, a rigorous derivation 
of the ad hoc expressions used by previous authors. The third model in this set can 
be applied to arbitrary flows, and it is stable because the C it predicts is always pos- 
itive. Finally, a model involving the subgrid-scale kinetic energy is presented which 
attempts to model back-scatter. This last model has some desirable theoretical fea- 
tures. However, even though it gives results in LES that are qualitatively correct, 
it is outperformed by the simpler constrained variational models. It is suggested 
that one of the constrained variational models should be used for actual LES while 
theoretical investigation of the kinetic energy approach should be continued in an 
effort to improve its predictive power and to understand more about back-scatter. 

2.1 A variational formulation 

Equation (4) may be written as E t j(x) = 0, where 

Eij(x) = Lij - ±L kk 6 t} - a,jC + %~C. (5) 

The residual Eij(x) at any given point depends on the value of the function C at 
neighboring points in the field. One cannot, therefore, minimize the sum of the 
squares of the residuals E X} E t j locally (as in Lilly, 1992) since reducing the value 
of EijE x j at one point changes its values at neighboring points. However, the 
method of least squares has a natural generalization to the non-local case. The 
function C that “best satisfies” the system of integral equations (4) is the one that 
minimizes 

T[C] = j £,,(x)£,,(x)<*x. (6) 

E[C] is a functional of C, and the integral extends over the entire domain. To find 
the Euler-Lagrange equation for this minimization problem, we set the variation of 
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T to zero: 

6E = 2 / Eij(x)6Eij(x)dx = 0. (7) 

Using the definition of Eij, we get 

j (-aijEijSC + E xi %SC) dx = 0 ( 8 ) 

which may be rearranged as 

J ctijEij + Pa J £ij(y)G(y,x)dy^ 6C(x)dx = 0. (9) 

Thus, the Euler- Lagrange equation is 

-otijEij + fcj J Eij( y)G(y,x)dy = 0 (10) 

which may be rewritten in terms of C as 

/(x) = C(x) - J fC(x,y)C(y)dy (11) 

where 

/(x) = — TT — 7~\ f«y( x ) £ <i( x ) - &;( x ) / Lij(y) G (y,x-W , 
a kl {x)a kl (x) [ J J 

r( V _ AC^(x,y) + ^(y,x) - K$(x,y) 

MX,yJ " « w (x)a«(x) 

and 

Ka(*, y) = a tJ (x)0ij(y)G(x,y), 

Ks(x, y) = Ai(x)A->(y) / &G(»,x)G(.,y). 

Equation (11) is readily recognized as Fredholm’s integral equation of the second 
kind. 

2.2 . The constrained variational problem 

In this section, we address the stability problem created by the negative eddy- 
viscosity by requiring that in addition to minimizing the functional (6), C satisfy 
some constraints designed to ensure the stability of the model. The choice of such 
constraints is clearly not unique. It is shown that the local least squares method 
(Lilly 1992) coupled with the volume averaging prescription (Germano et ai 1991) 
can actually be derived as a rigorous consequence of such a constrained variational 
problem for flows with at least one homogeneous direction. The method is then 
extended to general inhomogeneous flows. 
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2.2.1 Homogeneous turbulence 

In the case of homogeneous turbulence, it is natural to assume that C can depend 
only on time. Let us, therefore, impose this as a constraint in the problem of 
minimizing the functional (6). The functional F[C] then reduces to the function 

HC) = (CijCi,) - 2 (C, } m i} )C + (mijmjj)C 2 ( 12 ) 

where Cij ~ L{j — (l/3)L kk 8ij is the traceless part of rriij = — f3{j and ( ) 

denotes integral over the volume. The value of C that minimizes the function T(C) 
is easily found to be 

C = (L. } m,j) ( 13 ) 

(m kl m k i) 

where the isotropic part of Cij has vanished on contracting with the traceless tensor 
rriij. Equation (13) is precisely the result of Germano et al. and Lilly. 

2.2.2 Flows with at least one homogeneous direction 

As an example, we consider a channel flow with the y-axis along the cross-channel 
direction and periodic boundary conditions in the x and 2 directions. Since the flow 
is homogeneous in the x-z plane, we impose the constraint that C can depend only 
on time and the y co-ordinate. It is necessary to assume (as did Germano et al.) 
that the filtering kernel £?(x,y) is defined so as to be independent of the cross- 
channel direction, y. Therefore C may be taken out of the filtering operation and 
the functional (6) reduces to 

F[C] = J dy{(Cij - m,jC)(C,j - mijC)) xz (14) 

where ( ) xz denotes integral over the x-z plane and i = 1, 2, and 3 represents the 
x, y, and z directions, respectively. The condition for an extremal of the functional 
(14) may be written as 


f>T = 2 J dySC(y){m l} mijC - mij£,j) xz = 0 

which implies 

{rriij Cij rnijTrixjC) xz ~ 0 

and since C is independent of x and z and rriij is without trace, 

C — xz 

(mjt/mjt/) 

xz 


(15) 

(16) 


(17) 


This is the same expression as that of Germano et al. and Lilly for flows homoge- 
neous in the x-z plane. 
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2.2.3 Inhomogeneous flows 

In this section, we will adopt the point of view that perhaps the eddy-viscosity 
describes only a mean flow of energy from large to small scales, and back-scatter 
needs to be modeled separately as a stochastic forcing (Chasnov 1990, Leith 1990, 
Mason et al. 1992). We shall, therefore, insist on the eddy-viscosity always being 
positive and for the time being disregard back-scatter. Accordingly, in the problem 
of minimizing the functional (6), we impose the constraint 

C> 0. (18) 

It is convenient to write the variational problem in terms of a new variable £ such 
that C = £ 2 . Then the constraint (18) is equivalent to the condition that £ be real. 
In terms of the new variable £, equation (9) becomes 

j (- ctijEij + Pij j E tj (y)G{y,x)dyj £(x)^(x)dx = 0, (19) 

which gives for the Euler-Lagrange equation 

(-a i} Eij + 0ij j E i j(y)G(y,x)dy S j f(x) - 0. (20) 

Therefore, at any point x, either £(x) = 0 or the first factor in (20) vanishes. That 
is, at some points of the field C(x) = 0 and at the remaining points 

C(x) = G[C(x)} 


where 

£[C(x)] = f{x) + J IC(x,y)C(y)dy 

with f(x) and £(x, y) as defined in section 2.1. Note, however, we do not know 
in advance on which part of the domain C vanishes; this information is part of the 
solution of the variational problem. Therefore, if a C can be found such that 


rfic) = I Gl^(x)], if a[C'(x)J > 0; 
' ' 1 0, otherwise 


( 21 ) 


then it is a nontrivial solution of the Euler-Lagrange equation (20). Equation (21) 
may be written concisely as 


C(x) 


f(x) + J £(x,y)C(y)dy 


(22) 


where the operation denoted by the suffix *+’ is defined as + |^|) for 

any real number x . It is clear that a solution of (22) satisfies the Euler-Lagrange 
equation (20), but it is not obvious whether this solution is unique (we exclude the 
trivial solution C(x) = 0). Equation (22) is a nonlinear integral equation, and no 
rigorous results regarding the existence or uniqueness of its solutions are known to 
the authors. Nevertheless, we will assume that it has a unique solution in all cases 
of interest. Numerical experiments so far have given us no reason to question this 
assumption. 
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2.3. A model with back-scatter 

The instability associated with the negative eddy-viscosity may be understood 
in the following way. The Smagorinsky eddy-viscosity model does not contain any 
information on the total amount of energy in the subgrid scales. Therefore, if 
the coefficient C becomes negative in any part of the domain, the model tends 
to remove more energy from the subgrid scales than is actually available, and the 
reverse transfer of energy does not saturate when the store of subgrid- scale energy 
is depleted. However, in a physical system, if all the energy available in the subgrid 
scales is removed, the Reynolds stress will go to zero, thus quenching the reverse 
flow of energy. Clearly, a more elaborate model that keeps track of the subgrid-scale 
kinetic energy is required. Such a model is described in this section. (The possibility 
of treating the dynamic model in conjunction with an equation for turbulent kinetic 
energy was considered by Wong (1992) in a different context.) 

From dimensional analysis, the turbulent viscosity is the product of a velocity 
and a length-scale. We will take the square root of the subgrid-scale kinetic energy 
for the velocity scale and the grid spacing as the length scale. Thus, 



Tij - l 5ijT kk = — 2C'Afc 1 / 2 S, J 

(23) 

and 

Tij - ^ SijTtk = -2CAK l ' 2 %j 

(24) 

where 

k = ±( -«,«.) = 

1 

2 r "’ 

(25) 


A - -{uiUi - u,ui) - 

-T 

2 1 ' 1 ' 

(26) 

On taking the trace of (1), we have 




h — k + 2^”* 


(27) 


Since the average of the square of any quantity is never less than the square of its 
average, it follows that L lt is non-negative provided the filtering operation involves 
a non-negative weight G(x,y). Therefore, K is never less than fc, a result that 
might be anticipated since there are more modes below the test level cut-off than 
below the grid level. Substituting (23) and (24) in (1) and solving the corresponding 

variational problem, we get (11) with = — 2A I\ l ^ 2 Sij and j3 tJ = — 2A fc 1//f2 5, ; to 
determine C(x). 

To complete the model, it remains to give a method for determining k. For this 
we will use the well known model of the transport equation for k (e.g. Speziale 
1991) 


P/2 


dtk + Ujdjk = -TijSij -C.-— + dj(DAk l/2 d } k) + Re~ l d }J k 


A 


(28) 
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with the grid spacing A taken as the length-scale appropriate for the subgrid-scale 
eddies. Here C, and D are non- negative dimensionless parameters, and Re is a 
Reynolds number based on molecular viscosity. The coefficients C* and D can be 
determined dynamically. For this purpose, one writes down a model equation for 
K which is identical in form to (28) with test-level quantities replacing grid-level 
quantities. One then requires that K and k obtained by solving the correspond- 
ing evolution equations be consistent with (27). This gives the following integral 
equations for determining C, and D: 


C.(x) = 


/.(x) + j AC,(x,y)C.(y)dy 


(29) 


and 


D(x) = 


/d(x) + J £o(x,y )D(y)dy . 


The derivation and notation are explained in the appendix. 


(30) 


2.3.1 Stability 

It will be shown that the model described above is globally stable, that is, the 
total energy in the large-eddy field remains bounded in the absence of external 
forces and with boundary conditions consistent with no influx of energy from the 
boundaries. Using the continuity and momentum equations for the large-eddy fields 
and the sub-grid kinetic energy equation (28), we derive 


J t Q J UiU,dV + I kdV^j =- j ^k*' 2 dV - Re~ l J (djU t )(d } Ui)dV (31) 

where the integral is over the region occupied by the fluid. Boundary conditions are 
assumed to be such that there is no net flux of energy from the boundaries of the 
domain so that the surface terms vanish. Note that the terms in TijS,j which appear 
as a source term for k and a sink for the resolved scales (if C > 0 and vice versa 
when C < 0) have cancelled out in equation (31), and we are left with the result 
that in the absence of externally imposed forces and nontrivial boundary conditions, 
the total energy in the large and small scales taken together must monotonically 
decrease as a result of molecular viscosity. Using the notation 


and 


E(t) = j uiUidV 

(32) 

II 

(33) 


we have by (31), E(t) + e(<) < £(0) + e(0) and since e(<) > 0 (see next section), 
E(t) < 22(0) + e(0). Thus, the energy in the large-eddy field cannot diverge even 
though the eddy- viscosity is allowed to be negative. 
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2.3.2 Realizability 

It is necessary to demonstrate that the k computed using (28) has the following 
property; k(x,t) > 0 at all points x at all times t if fc(x, 0) > 0. This condition 
is required because it is clear from its definition that k cannot be negative, and, 
indeed, the model cannot be implemented unless the non-negativeness of k can be 
guaranteed. This condition is part of and included in a more general condition of 
‘realizability’ (Schumann 1977, Lumley 1978) required of subgrid-scale models to 
be discussed later. It must also be pointed out that La is an intrinsically positive 
quantity only if the filtering kernel G(x, y) > 0. The most commonly used filters 
in physical space such as the ‘tophat’ filter and the Gaussian filter do meet this 
requirement while the Fourier cut-off filter does not. Therefore, the Fourier cut-off 
filter may not be appropriate in this context. 

Suppose that initially (£ = 0), k > 0 at all points. Let t = to be the earliest time 
for which k becomes zero at some point x = Xo in the domain. It will be shown that 
dtk(xo,to) > 0 which ensures that k can never decrease below zero. Integration of 
equation (28) over an infinitesimal sphere of radius e centered around x 0 gives after 
dividing by e 3 

^ J u n kd(T + CAfc 1/,2 |5| 2 - ^ k 3/ 2 + J (34) 

where v = Re~ l + DA\/k and dcr is an infinitesimal element of area on the surface 
of the sphere with h as the outward normal. Since k first becomes zero at x = Xo, 
Xo is a local minimum. Therefore, k = Vfc = 0 at l Xo\ and hence k ~ e 2 and 
Vfc ~ e inside the sphere. Therefore, every term on the right side of (34) is of order 
e or higher except for the last term which is of order one. Thus, on taking the limit 
e — > 0 in equation (34), we have 

dk 1 [ dk J 

— = lim -r / v—da. (35) 

dt o e 3 J dn 

Since k is a minimum at the point Xo, the right-hand side is positive. Therefore, k 
can never decrease below zero. Note that we have assumed that C remains finite 
as k — > 0. Indeed, (27) implies remains finite as k (and hence fl XJ ) goes to zero. 
Thus, in this limit, (11) reduces to 

_ Q, J (x)I,;(x) 

^ ; a k ,{x)a kl (x) 

which is finite. Also, in this proof we assumed that the second derivatives of k 
at xo are not all zero. The proof, however, can be easily extended to remove this 
restriction. 

The requirement that k be non-negative is contained in a more general set of 
properties of the tensor r t j. They are called realizability conditions and may be 
stated in several equivalent forms (Schumann 1977). Since the Reynolds stress T t j 


12 


5. Ghosal, T. S. Lund & P. Moin 


is a real symmetric tensor, it can be diagonalized where the diagonal elements r a , 
T fj and r 7 are real. The realizability conditions can be stated as 

r ai T fi ,T^> 0. (36) 


It will be noted that (36) implies 

k = ^ th = ^(t„ + r 0 + r 7 ) > 0. (37) 

Positivity of the turbulent kinetic energy is, therefore, a consequence of the more 
general conditions (36). 

The modeled Reynolds stress (23) is diagonal in a co-ordinate system aligned to 
the principal axes of the rate of strain tensor, and the diagonal elements axe 

Ti = -2CAk l ^ Si + \k (38) 

where s; (t = a, /?, 7) are eigenvalues of the rate of strain. The realizability 
conditions (36) are satisfied if 

j fcl/2 Jfcl/2 

3A|s 7 | " “ 3As„ ' 


at each point of the field. In writing (39), the eigenvalues of the strain rate tensor 
have been arranged so that s Q > sp > s 7 . The incompressibiUty condition implies 
•St, + -S/j + s 7 = 0 and, therefore, s a > 0, s 7 < 0 and s 0 may be of either sign. 
Since C is obtained by solving the integral equation (11), it is difficult (perhaps 
impossible) to prove any general mathematical result on whether the realizability 
condition (39) is satisfied. Nevertheless, we offer the following estimates. 

A reasonable estimate for k when the turbulence is locally in equilibrium is k as 
A 2 |5| 2 . (This gives Smagorinsky’s formula on substitution in (23).) With this 
estimate for k, (39) may be written as C min < C < C max where 


V 2 + 4 + s 7 ^ v/2 

X |s 7 | “ 3 


(40) 


and 

c _ 

^ max — 0 ~ o ' V"**/ 

O o 

It is found from numerical experiments on LES of freely decaying homogeneous 
turbulence (section 2.5) that the distribution of values in the C field at any instant 
of time is approximately Gaussian with mean about 0.025 and standard deviation 
about 0.15. The precise values depend on the details of the simulation, but these 
numbers axe representative. For such a field, the number of points falling inside the 
range ^ — ^>nai exceeds 99.9 percent. 


y/*l + 4 + s * n /2 


(A 1 ^ 



A local dynamic model for LES 


13 


2.3.3 Galilean invariance 

The Navier- Stokes equations are invariant with respect to the transformations 

x ■ = xi — Uit (42) 

t = t (43) 

= Ui - Ui (44) 

where Ui is independent of space or time. It will be shown that the model described 
in this section as well as the constrained variational formulation lead to Galilean 
invariant equations for the large-eddy field. 

Substituting (44) in the definition of the Leonard term, we derive L- = X,j. 
Since Ui is constant, clearly the rate of strain tensors are invariant. It is shown 
in the appendix that equation (28) for the subgrid-scale kinetic energy is Galilean 
invariant. Therefore, k = fc, and on using equation (27), K = K. Therefore, each 
term in the integral equation for C in both the constrained variational formulation 
and the subgrid-scale kinetic energy formulation are invariant, so that C =C. The 
invariance of the model now follows on using -~r = and = jfa. 

2.3.4 Behavior near solid walls 

Consider a point near the wall with x \ , X2, and #3 in the streamwise, wall-normal, 
and spanwise directions, respectively. Then near the wall, u\ y u 1, U3, U3 ~ x 2 > and 
hence, by the continuity equation, U2, U 2 ~ (x 2 ) 2 * The near wall behavior of 
the subgrid-scale kinetic energy k involves some subtle issues. From its definition 
k = UiUi — UiTii , we must have k ~ (X2) 2 . In order to obtain such a behavior from 
(28), one needs to impose the boundary conditions that both k and dk/dx 2 vanish 
at the walls. However, since equation (28) is only second order in space, we cannot 
impose both these conditions. Thus, we are forced to choose only k = 0 at the wall 
and this, in general, will give a solution with the asymptotic behavior k ~ x 2. One 
possible remedy is to consider a two equation model (such as a k-e model) in place 
of equation (28). This gives a system that is fourth order in spatial derivatives and 
can, therefore, support the additional boundary condition (Durbin 1990. Also see 
the report by Cabot in this volume). In this report, however, we restrict ourselves 
to the simpler one equation model (28), and, therefore, the k obtained by solving 
(28) will in general have the behavior k ~ X 2 . Nevertheless, we will show that with 
the model coefficients computed dynamically, the eddy-viscosity is proportional to 
(X2) 3 and the molecular diffusion of kinetic energy balances the viscous dissipation 
near the wall independent of whether k ~ X 2 or k ~ (x 2 ) 2 . To stress this generality, 
we will write k ~ (2:2 ) 2m where m — 1/2 or 1 . 

The strain rate is dominated by the component S12 and S12 which are finite at 
the wall. The trace of the Leonard term La ~ (^2) 2 , hence by (27) K ~ (x2) 2m . 
Thus <*12, 0\2 ~ (^2 ) m are the only surviving terms of and fl t j near the wall. 
With these estimates, AC(x,y) ^ 1 so that 

C(x)w/(x)«^ ~(x 2 ) 3 - m . 

<*12 
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Therefore, the eddy-viscosity v t = CAy/Jc ~ (x 2 ) 3 , the well-known “y-cubed behav- 
ior at the wall”. 

We now consider the near wall behavior of the kinetic energy equation (28). Using 
the estimates given in the previous paragraph, ACp ~ 1, so that D ~ / D — > 0 at 
the wall. On using the expressions for AC* and /» given in the appendix, it is again 
easily verified that near the wall AC* ~ 1, and hence 

C.Jfe 3 / 2 ~ /.ik 3/2 ~ \ARe~ l djjLu\ 


which is finite. Hence, 


C,k 3 / 2 /A 

\Re~'djjk\ 


(45) 


independent of A and Re. Thus, all terms in the kinetic energy equation go to zero 
at the wall except the viscous dissipation and diffusion terms. These are finite and 
are of the same order at all Reynolds number and at any grid-resolution. This is 
the correct near wall behavior of the turbulent kinetic energy equation (Mansour, 
Kim and Moin 1988). 


2-4 Numerical methods 

The simplest iteration scheme for solving the integral equation (11) is 


C n+1 (x) = C„(x) + n 


/(x) + J AC(x,y)C„(y)dy - C n (x) 


(46) 


where ft is a relaxation factor that is selected to optimize convergence and n is an 
iteration index. To solve equation (22), the term 


/(x) + j AC(x,y)C'„(y)dy 

in (46) needs to be replaced by its positive part. Convergence can be achieved 
provided |/i| < fio where /i 0 is a number typically of order 0.1 for the present 
simulation. When the C(x) from the previous time-step was used to start the 
iteration, convergence at the level of one percent residual error took about twenty 
iterations. However, if a random velocity field is used in (46), values of n as small 
as 10~ 3 are required for stability, which greatly increases the number of iterations 
required for convergence. It is, therefore, prudent to use a better scheme that 
converges faster and is less dependent on the nature of the given LES field. Such a 
scheme can be devised using preconditioning methods. 

The integral equation (11) may be written in operator notation as 

(I - K)C = /. (47) 


On substituting K = E+(K-E) (where E is for the moment an arbitrary operator) 
in (47) we obtain 


C = (I - E)"7 + (I - E)-*(K - E)C. 


(48) 
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If C is replaced by C n +\ on the left-hand side and C n on the right-hand side we 
obtain the iteration scheme 


C„+i = (I - E )"7 + (I - E)-*(K - E)C„. (49) 


Equation (49) reduces to (46) if we choose E such that 


E ip = - — -ip 


(50) 


where t /> is an arbitrary function of position. It is well known that the speed of 
convergence of the scheme (49) depends on the eigenvalue spectrum of the operator 


B = (I — E)" 1 (K - E). 


(51) 


The smaller the spectral radius of B, the faster is the convergence. An efficient 
scheme is therefore obtained by choosing E such that ‘E « K’, and yet (I — E) can 
be readily inverted. One possible choice is 

EV’(x) = /iA 3 /C(x,x)V’(x) (52) 


where /i is a positive parameter. The motivation for (52) is the following; 

Kip = J fC(x,y)ip(y)dy = n(x)A* IC(x,x)ij>(x) (53) 

where fi(x) is expected to be between zero and one (corresponding to the limiting 
cases where the integrand is a delta function centered on x and a constant respec- 
tively). Equation (53) is exact. Equation (52) is the approximation to K obtained 
on ignoring the position dependence of /x. On substituting (52) in (49), we obtain 
after some algebra 


C„+i(x) = Y— /(x) + J £{x,y)C n (y)dy - ng(x)C„(x) 


(54) 


where 

y(x) = 


A 3 


a mn (x)a mn (x) 


2®i>(x)Aj(x)G(x,x) — 0ki(x.)0ki(x) J [G(z,x)] 2 dz 


(55) 

When G(x,y) is a tophat filter (that is, G(x,y) is constant inside the cube of edge 
A centered at x and zero outside), (55) reduces to 


fftophat(x) = 2a »i( X )/ ? »>( X ) ~ Ah(x)/Mx) 


Q'T7in(x)Q?frcn(x) 


(56) 
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The iteration scheme (54) can be modified to solve (22) on simply replacing 

/(x) + j >C(x,y)C n (y)dy 

by its positive part. The preconditioning scheme presented here is analogous to the 
point Jacobi scheme of inverting a matrix suitably generalized to the continuous 
case. Numerical tests indicate that (54) is considerably more robust than (46) and 
converges much faster. If ji is chosen between about 0.2 and 0.5, the scheme is found 
to be convergent for arbitrary velocity fields including fields of random numbers. 
With the C from the previous time step used as the initial guess, convergence at 
the level of 0.01 percent residual error never required more than three iterations for 
the LES of homogeneous isotropic turbulence described in the next section. 

2.5 Results 

In this section, results of LES of homogeneous isotropic freely decaying turbulence 
are presented using the constrained variational approach of section 2.2.3 (henceforth 
referred to as ieq-f ) and the kinetic energy approach of section 2.3. The results 
are compared to the results using the volume averaged (section 2.2.1) approach 
and to the grid-generated wind tunnel turbulence experiments of Comte-Bellot and 
Corrsin (1971). The Reynolds number based on the Taylor microscale and rms 
fluctuating velocity is in the range 40-70 in the experiment. The ratio of the integral 
to Kolmogorov scale is about 80. Therefore, even if the computational box is to be 
only three integral scales on a side, a DNS of the experiment will require at least 
(512) 3 grid points. 

A pseudo-spectral code due to Rogallo (1981) for homogeneous turbulence is 
modified and used with a (32) 3 grid. The filtering kernel G(x,y) is taken as a 
Gaussian function of width A = 2A, where A is the grid spacing. The results 
of the simulation are expected^ to be quite insensitive to the exact shape of the 
filtering kernel or of the ratio A/A (Germano et al 1991). However, the filtering 
kernel is restricted to be non-negative in the kinetic energy formulation because of 
considerations of realizability (section 2.3.2). 

Figure 1 shows the resolved energy decay as a function of time computed using 
the volume averaged, ieq+, and kinetic energy versions of the variational method 
together with the experimental data of Comte-Bellot and Corrsin (1971). All vari- 
ables are in non-dimensional units. The energy has been scaled with the mean 
square fluctuation velocity, (u /2 ) at the first measuring station. Time is in units of 
M/U where M is the grid size in the experiment and U is the free stream velocity. 
To facilitate comparison between computation and experiment, the experimental 
data has been filtered to remove spatial scales below the resolution of the compu- 
tational grid. Further, measured data on the energy decay as a function of distance 
from the grid has been converted into energy decay as a function of time using the 
measured convective velocity in the flow and Taylor’s hypothesis of 'frozen turbu- 
lence’. Both the volume averaged and constrained versions give results that are in 
good agreement with the experimental data. The result obtained using the kinetic 
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Time, t/MU~ l 

FIGURE 1. Decay of the resolved energy with time computed using the volume 

averaged ( — ), ieq+ ( * • • ) and kinetic energy ( ) versions of the variational 

method. 


energy method is seen to be not as good as the other two methods though it has 
the correct qualitative behavior. 

Figure 2 (A), (B), and (C) show the energy spectra at the initial time and two 
subsequent times computed using the volume averaged, ieq+, and kinetic energy 
versions of the variational method, respectively, together with the experimental 
points. Distances have been scaled by L/2n where L is length of the computational 
box which is about ten times the size of the grid generating the turbulence in the 
experiment. The initial velocity field is chosen to match the initial energy spectrum; 
therefore, only the last two curves on each figure are of significance. 

It is rather surprising that the kinetic energy method which apparently is based 
on more detailed physical ideas is outperformed by the other two relatively simple 
models els far as practical LES is concerned. In an effort to understand the reason 
for this, the subgrid-scale kinetic energy computed using (28) is plotted against 
the subgrid-scale kinetic energy derived from the experimental data. The result is 
shown in figure 3. The agreement with the experiment is seen to be remarkably 
good. Thus, we conclude that the relatively poor performance of the kinetic energy 
version of the model is not due to errors in prediction of the subgrid-scale kinetic 
energy. Perhaps the attempt to represent back-scatter by a negative eddy- viscosity 
is the source of the problem. These issues are still under investigation by the 
authors. 

A typical C(x) field computed using the constrained variational approach is shown 
in figure 4 over a horizontal cross section of the computational box. The field is 
seen to be highly variable with C fluctuating by as much as an order of magnitude 
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over a few grid spacings. Similar plots of C obtained by solving the (unconstrained) 
variational problem show a similar qualitative structure. 


3. Future plans 

The methods for solving the integral equations presented here represent only a 
first step. A systematic study of the numerical methods available for solving such 
integral equations must be undertaken. The numerical codes used to conduct the 
tests presented here have not yet been optimized for performance. This must be 
done before computations of flows with complex geometry are undertaken. The 
application of the method in computations with nonuniform grids involve some 
subtle issues that are being analyzed. It is expected that the methods presented 
here will be tested in progressively more challenging flows and the outcome of these 
tests will be used to guide future improvements. 

The kinetic energy version of the variational method is applicable to inhomo- 
geneous flows and provides an interesting model of back-scatter. However, the 
performance of this model in practical LES was not as satisfactory as that of the 
constrained variational methods. An attempt to understand the reasons for this is 
likely to produce not only a better model for LES but also a better understanding 
of back-scatter. Research in this direction is continuing. 

In conclusion, it must be stressed that the dynamic localization method is a 
very general idea and is not restricted in scope to Smagorinsky’s model or even 
to algebraic closure models. Extension of the ideas presented here to improved 
subgrid-scale models and higher order closures are interesting possibilities. 
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Appendix. The dynamic determination of C, and D. 

The evolution of subgrid- scale kinetic energy is modeled by the equation 

_ 1.3/2 

dtk + Ujdjk = -TijSij - C. — + d](DAk l / 2 djk) + {Re)~ l djjk (A.l) 

where C* and D are non-negative parameters. By analogy one can write down the 
corresponding evolution equation for the subtest-scale kinetic energy 

^ t"3/2 ^ 

d t K + %djK = -TijSij - C.^j- + dj(DAK^ 2 djK) + {Re)~'djjK. (A. 2) 

Equations (A. 1) and (A. 2) must be consistent with the relation 

K = k + ±Lu (A3) 

(equation (27) in the text). This fact can be exploited to determine C* and D as 
shown below. 

On substituting equation (A. 3) in (A. 2) one obtains after some algebraic manip- 
ulations ^ ^ 

d t k + Ujdjk s= -e + djfj + Re~ l djjk (A.4) 

where 

n r A ' 3 / 2 1 1 ^ 

e = TijS tJ + ^ -Re-'djjLii + -(d t Lu + UjdjLu) (A. 5) 

A ^ * 

and 

fj = DAK 1/2 djK. (A. 6) 

On applying the 1 A ’ operator to both sides of equation (.4.1) we obtain 

d t k -f %jdjk = -E + djFj 4- Re~ l djjk (.4.7) 


where 


and 


E = TijSij + (C,fc 3 / 2 /A) 


Fj = kuj — kuj + DAk l t 2 djk. 
Equations (A.4) and (A. 7) are consistent only if 


e = E 


(A.8) 

(A.9) 


(A.10) 


and 


fi = ^ 


(A.ll) 
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It may appear that the term | Re 'd }J L a has been arbitrarily lumped with e 
while it might as well have been converted to a divergence term and made part of 
fj . There is, however, no ambiguity as is evident from the following considerations. 
The coefficient D is a model for pressure diffusion and third order velocity correla- 
tion terms which do not depend explicitly on the Reynolds number. On the other 
hand, (7* models viscous dissipation which contains the Reynolds number explicitly. 
Thus, any term that contains the Reynolds number must appear in equation (A 10) 
which determines C* rather than in equation (A. 11) which determines D. A similar 
objection can be raised for the term (1/2 )ujdjLu in equation (A.5) which could 
have been written instead as -(1/2 )ujLa on the right hand side of equation (A. 6). 
The choice in this case is dictated by the requirement that both C* and D must 
be Galilean invariant in order to be consistent with their physical interpretation. 
Since it is only the total time derivative D t = d t + Ujdj appearing on the right hand 
side of equation (A. 5) and not d t that is a Galilean invariant operator, the apparent 
ambiguity is removed. Finally, the purist might argue that the curl of an arbitrary 
vector can be added on the right hand side of equation (A. 11) since we only need 
to enforce djfj = djFj. This degree of freedom merely refers to the fact that the 
definition of flux is always ambiguous up to the curl of an arbitrary vector field. 
Since we assume that the flux is defined identically at both the test and grid levels, 
the arbitrary solenoidal vector on the right hand side of equation (A. 11) must, in 
fact, be zero. 

Equations (A. 10) and (A. 11) can be rewritten in a more transparent form: 


X = (A. 12) 

and 

Zj = XjD - YjD (A. 13) 

with the notation 

X = T >jSij ~~ TijSij — —dtLn — -UjdjLu + - Re djjLn 



Zj 



= kuj — kuj 


Xj = AK 1/2 djK 

Yj = Ak'^djk. 


Clearly x, 4>, Zj, Xj, and Yj are known fields at any given time step. 

Since equation (A. 12) is a single integral equation for C,. one might be tempted 
to solve it directly without first going through the variational formulation. However, 
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since viscous dissipation can only remove energy from the subgrid scales, C m must 
be non-negative whereas the C* obtained by solving (A 12) can have either sign. 
One must, therefore, try to “best satisfy” (>1.12) subject to the constraint C* > 0. 
That is, C* must be chosen so as to minimize the functional 

j( x -<t>c.+w;) 2 dy 

with the constraint C* > 0. Similarly, D is obtained on minimizing 

J (Zj - XjD + Y~D)(Zj - XjD + Y)D)dy 

subject to the constraint D > 0. The solutions to these variational problems can be 
immediately written down in analogy to the solution (22) in the text if one notices 
the similarity in structure between equations (A, 12), (A. 13), and (4) in the text: 


where 


and 


C.(x)= /*(x) + J IC*(x,y)C*(y)dy 


/*( x ) = ^2^) «K x )x( x ) “ V’(x) J x(y)G(y,x)dy 

r /„ V N ^( x 1 y) + ^(y.x) -E$( x ,y) 

* ( ,y) 4>(*)<Hx) 

^( x .y) = «K X M y)G(x,y), 

£$( x » y) = 0( x Wy) J dzG(z,x)G(z,y). 


(A14) 


D(x) = 


/d{x) + J X D (x,y)D(y)dy (A.15) 


where 


and 


/D(X) = Xj(x)Xj(x) [ X >W Z M ~ Yj(x) J Z j(y)G{y ,x)dy , 

r /_ v x _ ^5( x -y) + £jt(y> x ) - £s( x >y) 

^-D(x,y) — Xj(x)Xj(x) 

£%(*, y) = *>( x )*>(y)G(x,y), 

IC$(x,y) = Y i (x)Y i (y) J dzG{z,x)G( z,y). 
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It is readily verified that x> V’, Z>> X } , find Yj are all Galilean invariant. 
Thus C* and D obtained by solving equations (A. 14) and (A. 15) are also Galilean 
invariant which, in turn, implies that the subgrid-scale kinetic energy equation (A.l) 
itself is Galilean invariant. 
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